home *** CD-ROM | disk | FTP | other *** search
/ PC-X 1997 October / pcx14_9710.iso / swag / math.swg / 0115_Calculating a Factor using Gamma.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1996-08-30  |  2.2 KB  |  75 lines
  1. (*
  2.                 Daniel Doubrovkine - dblock@infomaniak.ch
  3.    part of the Expression Calculator 2.0 Multithread, destribute freely
  4.               http://www.infomaniak.ch/~dblock/express.htm
  5.   (ref. Wanner / Hairer - Analysis by it's History - University of Geneva)
  6.  
  7.   calculating a factor of any defined value inclusing negative and positive
  8.       non integers! using Euler's (another swiss guy) Gamma function:
  9.        n!=Gamma(n+1) and (n-1)!=Gamma(n)/(n-1) for negative values
  10.         Gamma(Alpha):=Integral zero->infinity of x^(Alpha-1)*E^-x
  11.  
  12.  in GammaIntegral function TOTAL PRECISION (not just approx!) is reached by
  13.    calculating an integral from 0 to 100 with a step if 0.01 from Gamma(4)
  14. *)
  15. function Gamma(alpha,step:extended):extended;
  16.  function GammaIntegral(alpha:extended):extended;
  17.   (*x^y:=e^(y*ln(x));*)
  18.   function Power(base, exponent: extended): extended;
  19.    begin
  20.       Power:=exp(exponent*ln(base));
  21.       end;
  22.   function IntegralStep(x: extended):extended;
  23.   begin
  24.      (*Gamma Integral Step...just to have less mess*)
  25.      IntegralStep:=power(x,alpha-1)/power(Exp(1),x);
  26.      end;
  27.  (*Gamma Integral*)
  28.  var
  29.    GammaTempIntegral:extended;
  30.    l: extended;
  31.  begin
  32.    l:=0;
  33.    GammaTempIntegral:=0;
  34.    while l<100 do begin
  35.          l:=l+Step;
  36.          GammaTempIntegral:=GammaTempIntegral+IntegralStep(l)*step;
  37.        end;
  38.        GammaIntegral:=GammaTempIntegral;
  39.        end;
  40. (*Gamma*)
  41. var
  42.   NewGamma: extended;
  43.   i: integer;
  44.   begin
  45.   if (alpha<=0) then begin
  46.     if (trunc(alpha)=alpha) then begin
  47.        writeln('factor results to infinite value');
  48.        halt;
  49.        end
  50.        else begin
  51.             NewGamma:=GammaIntegral(abs(1+frac(alpha))+1);
  52.             for i:=0 to trunc(abs(alpha))-1 do begin
  53.                 NewGamma:=NewGamma/(frac(alpha)-i);
  54.                 end;
  55.             Gamma:=NewGamma;
  56.             exit;
  57.             end
  58.   end  else begin
  59.      (*the following values return trailing decimals, this corrects it*)
  60.      if alpha=1 then Gamma:=1 else
  61.      if alpha=2 then Gamma:=1 else
  62.      if alpha=3 then Gamma:=2 else
  63.      Gamma:=GammaIntegral(alpha);
  64.      end
  65. end;
  66.  
  67. function Factor(n: extended):extended;
  68.  begin
  69.    Factor:=Gamma(n+1,0.01);
  70.    end;
  71.  
  72. begin
  73.      writeln(Factor(6):0);
  74. end.
  75.